##############################################################################
# Filename: Analysis_DemonetizationOpinion.R
# Purpose: Produce Table 4 and Figure 3 in main text and SI Figure A2
##############################################################################

source("Setup.R")

### In-text discussion of group means
data_op %>%
  group_by(Article, ProBJP) %>%
  summarise(
    count = n(),
    mean = mean(as.numeric(op_dem), na.rm = TRUE),
    sd = sd(as.numeric(op_dem), na.rm = TRUE)
  )


### Table 4: Effect of article read and level of BJP support on 
          #  opinion about demonetization
# Level of agreement with statement
mod_dem_main = lm(as.numeric(op_dem) ~ ProBJP*Article
                    + CollegeGrad
                    + NewsDaily + StrongInterestPolitics,
                    data = data_op)
# Certainty (neither agree nor disagree)
data_op$op_dem_neither = ifelse(data_op$op_dem=="Neither agree nor disagree",1,0)
mod_dem_neither = lm(op_dem_neither ~ Article_Alignment 
                             + CollegeGrad 
                             + NewsDaily + StrongInterestPolitics,
                             data = data_op)
# Print
stargazer(mod_dem_main, mod_dem_neither,
          align=TRUE, keep.stat="n", no.space=TRUE,
          covariate.labels=c("Pro-BJP", "Negative Article", "Positive Article",
                             "Counter-Attitudinal Information", "Pro-Attitudinal Information",
                             "College Graduate", "Daily News",  "Interested in Politics",
                             "Pro-BJP x Negative Article", "Pro-BJP x Positive Article")
)


### Figure 3: Relationships between level of BJP support and 
           #  demonetization opinion by article 
plot.data = data_op[,c("ResponseId","ProBJP","Article","op_dem")]
plot.data$op_dem = as.numeric(plot.data$op_dem)
plot.data = reshape2::melt(plot.data, id.vars = c("ResponseId","ProBJP","Article"))
plot.data$variable = NULL
plot.data$Article = factor(plot.data$Article, levels=c("DN","pl","DP"), 
                    labels=c("Negative Article","Control","Positive Article"))
plot.data = plot.data %>% 
  group_by(Article, ProBJP) %>%
  summarise(mean = mean(value,na.rm=TRUE),
            sd = sd(value,na.rm=TRUE),
            n = length(value),
            median = median(value,na.rm=TRUE))
plot.data$se = plot.data$sd / sqrt(plot.data$n)  
pd = position_dodge(.6)   
ggplot(plot.data, aes(x = ProBJP, y = mean, 
                      group = Article, color = Article, shape = Article)) +
  geom_point(size = 1, position = pd) +
  geom_errorbar(aes(ymin  = mean - se,
                    ymax  = mean + se),
                width = 0.1,
                size  = 0.4,
                position = pd) +
  scale_y_continuous(name="Mean",limits = c(1,5)) +
  scale_color_manual(values = c("darkblue","darkviolet","darkgreen"),
                     name = "") +
  scale_shape_manual(values = c(15,17,16),
                     name = "") +
  coord_fixed(ratio = .5) +
  theme_bw() + theme_minimal() + 
  theme(plot.title = element_text(family = "sans", size=13, hjust=0.5),
        axis.title = element_text(family = "sans", size=11),
        axis.text.x = element_text(family = "sans", size=10),
        legend.position= "bottom",
        plot.background = element_rect(fill="white",color="white"))
# ggsave("plot_int_OpDem.png", width = 5, height = 4)
rm(plot.data,pd)


### Figure A2: Partisan gap in issue opinions
data_controlgroup = data_op[which(data_op$Article=="pl"),]
opinion_vars = c("op_fake","op_ubi","op_res","op_gst","op_dem","op_kash","op_dev","op_min","op_imm","op_women")
bjp_est = numeric()
bjp_ci_low = numeric()
bjp_ci_upp = numeric()
bjp_pval = numeric()
for (i in 1:10){
  themodel = lm(as.numeric(data_controlgroup[,opinion_vars[i]]) ~ ProBJP + CollegeGrad + NewsDaily + StrongInterestPolitics + Age + Hindu + Male, data = data_controlgroup)
  bjp_est[i] = coefficients(themodel)[2]
  bjp_ci_low[i] = confint(themodel)[2,1]
  bjp_ci_upp[i] = confint(themodel)[2,2]
  bjp_pval[i] = summary(themodel)$coefficients[2,4] 
}
opinion_vars = c("Fake news","UBI","Reservations","GST","Demonetization","Kashmir","Development schemes","Minorities","Immigrants","Women in politics")
bjp_ests = data.frame(opinion_vars, bjp_est, bjp_ci_low, bjp_ci_upp, bjp_pval)
bjp_ests = bjp_ests[order(bjp_ests$bjp_est),]
bjp_ests$opinion_vars = factor(bjp_ests$opinion_vars, levels = bjp_ests$opinion_vars)
ggplot(bjp_ests, aes(x = opinion_vars, y = bjp_est)) +
  geom_hline(aes(yintercept=0), size = .5, color="gray") +
  geom_errorbar(aes(ymin=bjp_ci_low, ymax=bjp_ci_upp), size=.2, width=.5) +
  geom_point( size=2, color="black") +
  scale_color_manual(values=bjp_ests) + 
  labs(x = "", y = "") +
  coord_flip() +
  theme_bw() + theme_minimal() 


### In-text discussion of differential levels of opposition to demonetization  
data_controlgroup$oppose_dem = ifelse(data_controlgroup$op_dem=="Disagree"|data_controlgroup$op_dem=="Strongly disagree",1,0)
# Among the 285 strong BJP supporters, 216 (24%) oppose demonetization;
# among the 202 non-BJP supporters, 154 (76%) oppose demonetization
table(data_controlgroup$oppose_dem, data_controlgroup$ProBJP)